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"55 ■ ABSTRACT 
>v 

We present rigorous 3D EMF simulations of isolated features on photomasks using a newly developed finite- 
Qh element method. We report on the current status of the finite-element solver JCMsuite, incorporating higher- 
order edge elements, adaptive refinement methods, and fast solution algorithms. We demonstrate that rigorous 
, and accurate results on light scattering off isolated features can be achived at relatively low computational cost, 
' compared to the standard approach of simulations on large- pitch, periodic computational domains. 
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1. INTRODUCTION 

Rigorous electromagnetic field (EMF) simulations have become an important tool for mask design in low k\ 
, applications. As has been shown in previous works finite-element simulations are superior in terms of simulation 
accuracy, convergence rate and computation speed when compared to other presently used EMF simulation 
^ . methods. 1 ' 2 We address 3D simulation tasks occuring in microlithography by using the frequency-domain 
FEM solver JCMsuite. This solver has been successfully applied to a wide range of 3D electromagnetic field 
computations including microlithography, 1-4 left-handed metamaterials in the optical regime, 5 and photonic 
crystals. 6 The solver has also been used for pattern reconstruction in EUV scatterometry, 7 ' 8 and it has been 
benchmarked to other methods in typical DUV lithography projects. 1,2 

In this paper we report on the current status of the finite-element solver JCMsuite, and we present simulation 
results of light transition through isolated and periodic features on photomasks. 
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2. BACKGROUND 



The finite-element package JCMsuite allows to simulate a variety of electromagnetic problems. It incorporates a 
scattering solver (JCMharmony) , a propagating mode solver (JCMmode) and a resonance solver (JCMresonance) . 
The scattering, eigenmode and resonance problems can be formulated on ID, 2D and 3D computational domains. 
Admissible geometries can consist of periodic or isolated patterns, or a mixture of both. Further, solvers for 
problems posed on cylindrically symmetric geometries are implemented. 
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Figure 1. Schematics of a mask with an isolated line: cross-section in a x-y- plane. The line with a critical dimension of 
CD X is depicted in light grey. The pattern extends towards infinity in positive and negative x- and y-directions. 
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Figure 2. Schematics of an isolated line mask: cross-section in the z-z-plane. The material stack on a Quartz (Si02) 
substrate consists of two layers containing chromium. The computational domain is indicated. Please note that the 
computational domain is significantly smaller than the computational domain for a periodic line mask with high p x /CD x - 
ratio (see Fig. 4). 

In this paper we concentrate on light-scattering off 2D and 3D photomask patterns (lines and contact holes). 
The patterns of interest are isolated in the x— and y— directions and are enclosed by homogeneous substrate 
and superstrate (typically air) which are infinite in the —z-, resp. +z-direction. Cross sections through these 
patterns are schematically shown in Figures 1 and 2. We compare the results obtained from calculations on 
isolated computational domains to results from calculation on periodic computational domains with a large 
pitch. Cross sections through the periodic patterns are schematically shown in Figures 3 and 4. 

Light propagation in the investigated systems is governed by Maxwell's equations where vanishing densities of 
free charges and currents are assumed. 9 The dielectric coefficient e(x) and the permeability /x(x) of the considered 
photomasks are complex and in the case of periodic mask layouts periodic, e{x) = e (x + a), n(x) = (x + a). 
Here a is any elementary vector of the periodic lattice. For given primitive lattice vectors a\ and S2 a periodic 
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Figure 3. Schematics of a periodic line mask: cross-section in a x-y-plane. Lines with a critical dimension of CD X are 
depicted in light grey. The pattern is periodic in :r-direction with a pitch of p x . 
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Figure 4. Schematics of a periodic line mask: cross-section in the :r-z-plane. The material stack on a Quartz (Si02) 
substrate consists of two layers containing chromium. The computational domain is indicated. 



elementary cell f2 C K 3 can be denned as fi = {x <E M 2 | x = a\a\ + 0:202; < a\, a.<i < l} x [z su b, z sup \. In 
the isolated case the computational domain is choosen such that in the exterior domain only homogeneous or 
waveguide- like structures are present (see Fig. 2). 

A time-harmonic ansatz with frequency w and magnetic field H(x, t) = e~ lu]t H(x) leads to the following 
equations for H{x): 
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Figure 5. Schematics of the simulation flow of the FEM solver JCMsuite. 



layout setup 


isolated 


periodic 




set 1 


set 2 


set 3 


CD X [nm] 
d Cr,top I nm ] 
d Cr,bottom I nm ] 


260 

18 

55 


Px 

w x 


400 nm 


6 x A : 4 x A : 206 x A 


1 jjLva : 40 nm : 20 fim 


Ao 

e rair 

rCr,top 
e rCr, bottom 
e rSiO, 


193.0 nm , coherent illumination, a = 
1.0 

(1.965+ 1.20H) 2 
(1.477 + 1.762i) 2 
(1.564) 2 


N.A. 

Magnification 
Rcl. threshold 


1.0 
4 X 
0.69 



Table 1. Parameter settings for the isolated (data set 1) and periodic (data sets 2 and 3) line mask simulations: Mask 
geometry, material parameters, illumination, imaging parameters. The width of the computational window in the isolated 
case is w x = 400 nm (set 1), the pitch is varied in steps of 4 x Ao, resp. in steps of 40nm in the periodic case (set 2, 
resp. 3). The relative permittivity, e r , is the square of the complex index of refraction, e r = n 2 = (n r + ik) 2 . 



The wave equation and the divergence condition for the magnetic field: 

V x — !— V x H{x) - w 2 fi(x)H(x) = 0, xefl, (1) 

£{X) 

V • fj,(x)H(x) = 0, fed. (2) 

In the case of an isolated pattern: Transparent boundary conditions at all boundaries of the computational 
domain, Oil, where H m is the incident magnetic field (plane waves in this case), and n is the normal vector 
on dQ,: 

j-\7 x (H — H" 1 )) xn = DtNiH - H m ), x G dSl. (3) 

^£{X) 

The DtN operator (Dirichlet-to-Neumann) is realized with an adaptive PML method. 10,11 This is a 
generalized formulation of Sommerfeld's radiation condition. 

• In the case of a x — y-periodic pattern: Transparent boundary conditions at the boundaries to the substrate 
(at z su b) and superstrate (at z sup ), <9f2, according to Equation 3, and periodic boundary conditions for the 
transverse boundaries, <9f2, governed by Bloch's theorem 12 : 

H(x) = e iS ' s u(x) 7 u(x) = u(x + a), (4) 

where the Bloch wavevector k £ R 3 is defined by the incoming plane wave H m . 

Similar equations are found for the electric field E(x,t) — e~ lut E(x)\ these are treated accordingly. The 
finite-clement method solves Eqs. (1) - (4) in their weak form, i.e., in an integral representation. 

The finite-element methods consists of the following steps: 

• The computational domain is discretized with simple geometrical patches, JCMsuite uses linear (ID), 
triangular (2D) and tctrahcdral or prismatoidal (3D) patches. The use of prismatoidal patches is well 
suited for layered geometries, as in photomask simulations. 



• The function spaces in the integral representation of Maxwell's equations are discrctized using Ned- 
elec's edge elements, which are vectorial functions of polynomial order defined on the simple geometrical 
patches. 13 In the current implementation, JCMsuite uses polynomials of first (1 st ) to ninth (9 th ) order. In 
a nutshell, FEM can be explained as expanding the field corresponding to the exact solution of Equation (1) 
in the basis given by these elements. 

• This expansion leads to a large sparse matrix equation (algebraic problem). To solve the algebraic prob- 
lem on a standard workstation linear algebra decomposition techniques (LU-factorization, e.g., package 
PARDISO, 14 which was used in the simulations of Chapters 3, 4) are used. In cases with either large 
computational domains or high accuracy demands, also domain decomposition methods 15 are used and 
allow to handle problems with very large numbers of unknowns. 

For details on the weak formulation, the choice of Bloch-periodic functional spaces, the FEM discretization, 
and the implementation of the adaptive PML method in JCMsuite we refer to previous works. 10 ' 11 In future im- 
plementations performance will further be increased by using curvilinear elements, general domain-decomposition 
techniques and /ip-adaptive methods. 




Figure 6. (a) Dependence of aerial image CD on computational domain width w x , resp. pitch p x for isolated and periodic 
layouts according to Table 1. The single data point (+) corresponds to the isolated layout and solver definitions (data 
set 1), the data points on the dotted/solid line correspond to the periodic layout (data set 2/3). The strong variation of 
CD with pitch for the points on the dotted line is not observed when the CD is probed at specific pitch values only (solid 
line), (b) Detailed view of the same data as in (a), (data sets 1 and 3). 



3. RIGOROUS SIMULATIONS AND LARGE-PITCH PERIODIC-DOMAIN 
SIMULATIONS OF AN ISOLATED LINE 

We apply the light-scattering solver of the programme package JCMsuite in order to simulate the lithographic 
projection of an isolated line. 

The setup is as schematically shown in Figures 1 and 2, with geometrical, material, illumination and imaging 
parameters as in Table 1 (data set 1). We compare the obtained aerial image to results obtained with a setup 
of periodic lines of otherwise the same parameters, with a large ratio of pitch to linewidth (see Figures 3, 4 and 
Table 1, data sets 2, 3). 

Figure 5 shows the simulation flow of the finite-element software: 

• The geometry of the computational domain is described in a polygonal format corresponding to Fig. 2 with 
material attributions and statements about the computational domain boundaries (periodic and trans- 
parent boundaries in this case). The translational vectors of the periodic pattern (01,122) are identified 



automatically from the layout, optimized settings of the perfectly matched layers (PML) are found auto- 
matically with an adaptive method (adaptive PML, aPML). Further, parameters specifying a maximum 
patch size of the finite-element mesh and further meshing properties can be set. From these geometry 
parameters and mesh parameters the mesh is generated automatically. 

• Material parameters (complex permittivity and permeability tensors) can be specified as piecewise constant 
functions and/or (using dynamically loaded libraries) as functions of arbitrary spatial dependence. In the 
presented example, the piecewise constant, isotropic settings given in Table 1 are used. The relative 
permeability is pL r — 1 for all used materials. 

• Light sources can be defined as predefined functions (plane waves, Gaussian beams, point sources) or as 
arbitrary functions using dynamically loaded libraries. JCMsuite allows to generate solutions to several 
independent source terms in parallel, efficiently re-using the inverted system matrix. For the simulations in 
this chapter, without loss of generality, a coherent, polarized point source is used (wavevector k = (0, 0, k z )). 
For the simulations in Chapter 4, an unpolarized, incoherent source is used. 

• The main project definitions in this case are the accuracy settings (mesh refinement, PML refinement, 
finite-element degree) and the definitions of postprocesses. Upon execution the JCMsolve computes the 
full field distribution over the entire computational domain. Through internal or external post-processes, 
the quantities of interest can be derived from the field. E.g., the complex far field coefficients can be 
attained by Fourier integration / transformation and by the evaluation of the Raylcigh-Sommcrfeld diffrac- 
tion formula, 16 an aerial image is calculated, or the field distribution is exported to graphics format for 
visualisation and analysis. 

• Interfaces to scripting languages like Python or MATLAB can be used for performing automatic parameter 
scans and data analysis. In the presented example, the interface for the generation of the input files and for 
the execution of a loop over the geometrical parameters given in Tables 1, 2 has been realized in MATLAB. 
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Figure 7. Schematics of a mask with an isolated contact hole: cross-section in an x-y-pl&ne. The contact hole with a 
critical dimensions in x-, resp. {/-direction of CD X / CD y is depicted in light grey. The pattern extends towards infinity 
in positive and negative x- and y- direct ions. A cross-section through the computational domain is indicated. 

From the obtained aerial image intensity distribution we decude a mask CD of the printed feature. For the 
specific parameter setting as given in Table 1 and from a simulation on the isolated computational domain we 
obtain a well converged value of the CD of approximately 65 nm. When we instead simulate the light distribution 
using periodic computational domains with large pitches we obtain some dependence of the CD on the pitch of 
the chosen computational domain. As expected, with increasing ratio of pitch to linewidth, the CD resulting 
from the aerial image at focal position with a fixed threshold (0.69 of the maximum intensity in this case) is 
converging to a fixed value for pitches well above 5 microns. This behavior is depicted in Figure 6. Please note 
that due to across-pitch imaging effects the resulting CD variation is rather large for a quasi continuously varied 
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Figure 8. Schematics of a 2D-periodic array of contact holes: cross-section in a x-y-p\ane. Holes with critical dimensions 
of CD X , CD y are depicted in light grey. The pattern is periodic in x /y- direct ion with a pitch of p x /p y - 
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Table 2. Parameter settings for the isolated (data set 1) and periodic (data set 2) contact hole mask simulations: Mask 
geometry, material parameters, illumination, imaging parameters. The width of the computational window in the isolated 
case is w x = 400 nm x w y = 400 nm (set 1). The pitch is varied in steps of 40 nm in the periodic case (set 2). The source 
is an unpolarized, single-point source. 



pitch (see, e.g. 17 ), see the dotted line in Figure 6(a). A smoother dependence is observed when the pitch is 
probed at multiples of wavelength times imaging magnification, sec solid line in Figure 6 (a). However, also with 
this method, the waver scale CD error at a pitch-to-linewidth ratio of 10:1 is still of the order of 1 nm. Therefore, 
for accurate results on isolated features using a periodic model with large pitch-to-linewidth ratio, very large 
computational domains have to be chosen. 

This demonstrates the advantages using the non-periodic model according to Fig. 2 and Equation 3. Using 
this model, the computational domain can be chosen just slightly larger than the feature of interest. The 
simulated aerial image still coincides with the aerial image from the periodic model with largest investigated 
pitch. This makes it possible to accurately compute images of isolated features with low computational effort 
(see singular data point for the isolated case in Figure 6 a,b). 

4. 3D ISOLATED CONTACT HOLE - NEAR FIELD CONVERGENCE 



In this chapter we investigate the application of our method to the simulation of light scattering off a 3D isolated 
contact hole. We concentrate on the convergence of the near field and - as in Chapter 3 - on the comparison to 




Figure 9. 2D cross section through the 3D near field intensity distribution at the upper chromium layer. The colorbar 
on the right indicates the values for the field intensity / = I{x,y, zq) (arbitrary units). 
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Figure 10. ID cross section through the 3D near field intensity distribution. The red solid line corresponds to the 
isolated case (Table 2, data set 1). Black lines correspond to the periodic comparison case with different pitches as given 
in Table 2, data set 2. 

the large-pitch periodic case. Figure 7 shows a schematics of the geometrical layout: A contact hole with square 
shape of width CD X and length CD y (in mask scale dimensions) spans in the x-y-plane. A cross-section in the 
x-z-plane corresponds the lincmask case, as depicted in Figure 2. 
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Figure 11. Field intensity of the data from Figure 10 at x = 0. The data point for the simulation on an isolated domain 
is indicated by a red circle, the data points for simulations on periodic domains with varied pitches p x and p y = p x is 
indicated by crosses. 



Figure 9 shows a x — y-cross section through the obtained near field intensity distribution. The cross section is 
taken at constant z position of zq = —dcr,top/2, where z = corresponds to the upper edge of the upper chromium 
layer. The illumination in this case is an unpolarized, incoherent point-source at perpendicular incidence from the 
substrate. For the finite-element expansion and for the realisation of the PML boundary condition, fourth-order 
finite-elements have been used. The calculation has been performed on a standard 64-bit-CPU workstation with 
extended memory and several dual core processors (AMD Opteron). The computation time for this problem 
with about 100,000 unknowns is of the order of one minute. The convergence of the intensity 1(0, 0, zq) is plotted 
in Figure 12. Here, fourth order finite elements are used, and simulations with different numbers of unknowns 
in the FEM expansion are reached by different spatial discretizations of the computational domain. Figure 13 
shows the computational effort in minutes of total computation time in dependence of numbers of unknowns in 
the FEM expansion of the near field solution, N. The values give only an estimate of the computational costs 
as upon the computations, the workstation has been used for other tasks in parallel. Typical memory (RAM) 
usage for the FEM computation of the solutions corresponding to the data points of this Figure is between 1 
and 15 GigaByte. 

Figure 10 shows a ID cross section through the intensity I(x, y, z) of the near field solution at constant 
zq = —dcr,top/2 and y = 0. Additional to the result obtained on the isolated periodic domain, results obtained 
on periodic domains are plotted in this figure. Due to typical computational resources limitations of full 3D 
computations the pitch in this case can not be choosen very large. Therefore the variation of the obtained field 
distributions in the investigated parameter regime is still of significant magnitude. 

Figure 11 shows the intensity values corresponding to the data from Figure 10 at x = 0. As can be seen 
from the figure and as in the line mask case, a rather large computational domain has to be choosen when 
one tries to compute an accurate near field solution of an isolated contact hole on a periodic computational 
domain. The computational costs of the isolated and periodic computations as depicted in Figures 11 and 12 
are shown in Figure 13. As can be seen from the Figures, an accurate solution where the field intensity at a 
specified position is converged with an error of lower than 0.5% can be reached with about 100,000 unknowns 
on the isolated computational domain. In contrast, with a large-pitch periodic model, the pitches in x- and 
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Figure 12. Convergence of the near field intensity {I(ro) as function of the number of unknowns in the simulation). 

y-direction have to be chosen very large to ensure that the model-induced error is neglectible. This leads to a 
very high number of unknowns in the FEM problem and consequently to a highly increased computational effort. 
Therefore, and especially for more complex isolated problems under investigation, a large-pitch periodic model 
is not a feasible option. This also holds for other discretization methods than FEM using periodic models for 
isolated problems. For the demonstrated 3D test problem the advantage of using isolated computational domains 
with transparent boundary conditions on all boundaries is more than one order of magnitude in computation 
time. With the implementation of this method rigorous and accurate investigation of scattering from isolated 
features has become straight-forward. 

5. CONCLUSIONS 

We have performed rigorous 3D FEM simulations of light transition through isolated features on linemasks 
and contact hole photomasks using the finite-element programme package JCMsuite. We have investigated the 
convergence behavior of the solutions and we have shown that we achieve results at high numerical accuracy. We 
have compared the results to calculations on periodic, large-pitch domains and found a performance advantage 
(cpu time) of at least one order of magnitude for calculations on isolated computational domains. Our results 
show that rigorous 3D simulations of isolated features on masks or wafers or in other settings can well be handled 
at high accuracy and relatively low computational cost. 
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